Ab initio calculation of valley splitting in monolayer 5-doped phosphorus in silicon 
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The differences in energy between electronic bands due to valley splitting are of paramount im- 
portance in interpreting transport spectroscopy experiments on state-of-the-art quantum devices 
defined by scanning tunneling microscope lithography. We develop a plane-wave density functional 
theory description of these systems which is size-limited due to computational tractability. We then 
develop a less resource-intensive alternative via localized basis functions, retaining the physics of 
the plane-wave description, and extend this model beyond the capability of plane-wave methods to 
determine the ab initio valley splitting of well-isolated 5-layers. In obtaining agreement between 
plane-wave and delocalized methods, we show that the valley splitting has been overestimated in 
previous ab initio calculations by more than 50%. 
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I. INTRODUCTION 

Quantum devices in silicon have been the subject of 
concentrated interest, both experimental and theoreti- 
cal, in recent years. Efforts to make such devices have 
led to atomically precise fabrication methods which in- 
corporate phosphorus atoms in a single monolayer of a 
silicon crystalir— . These dopant atoms can be arranged 
into arrays^, or geometric patterns for wires^ and as- 
sociated tunnel junctions 7 -, gates, and quantum dots^ 
- all of which are necessary components of a function- 
ing device 1 ^. The patterns themselves define atomically 
abrupt regions of doped and undoped silicon. While sili- 
con, bulk-doped silicon, and the physics of the phospho- 
rus incorporation^ are well-understood, models of this 
quasi-two-dimensional phosphorus sheet are still in their 
initial stages. In particular, it is critical in many ap- 
plications to understand the effect of this confinement 
on the conduction band valley degeneracy, inherent in 
the band structure of silicon. For example, the degener- 
acy of the valleys has the potential to cause decoherence 
in a spin-based quantum compute r 12 ' 13 , and the degree 
of valley degeneracy lifting (valley splitting) defines the 
conduction properties of highly confined planar quantum 
dots^. 

The importance of understanding valley splitting in 
monolayer (5-doped Si:P structures has led to a number of 
theoretical works in recent years spanning several tech- 
niques, from pseudopotential theories via planar Wan- 
nier orbital (PWO) basest, density functional theory 
(DFT) via linear combination of atomic orbital (LCAO) 
base a 15 ' 16 , to tight-binding (TB) models^ 7 - - — , and effec- 
tive mass theories (EMT)22r— . We note that several of 
these papers are based upon the assumption that the 
effective masses of <5-doped P in Si remain unchanged 
from bulk-doped values j 22 ' 23 an assumption which has 
been challenged] 14 ' 17 Others assume doping over a multi- 
atomic plane ban d 17 ' 22 which no longer represents the 
state of the art in fabrication. We also note that Rcf. 
[l5l represents the first attempt to model these devices 



by considering explicitly doped 5-layers with DFT, using 
a relatively small localized basis set with the assump- 
tion that a basis set sufficient to describe bulk silicon 
would also adequately describe P-doped Si. There is cur- 
rently little agreement between the valley splitting values 
obtained using these methods, with predictions ranging 
between 5 to 270 meV, depending on the arrangement 
of dopant atoms within the 5-layer. Density functional 
theory has been shown to be a useful tool in predict- 
ing how quantum confinement and/or doping perturbs 
the bulk electronic structure in silicon- and diamond-like 
structures^ - — and it might be expected that the re- 
moval of the basis set assumption will lead to the best 
estimate of the valley splitting available. 

In this paper we determine a consistent value of the 
valley splitting in explicitly (5-dopcd structures by ob- 
taining convergence between distinct DFT approaches 
in terms of basis set and system sizes. We perform a 
comparison of DFT techniques, involving localized nu- 
merical atomic orbitals and delocalized plane- wave (PW) 
basis sets. Convergence of results with regard to the 
amount of Si "cladding" about the <5-doped plane is stud- 
ied. This corresponds to the normal criterion of super- 
cell size, where periodic boundary conditions may intro- 
duce artificial interactions between replicated dopants in 
neighboring cells. A benchmark is set via the delocal- 
ized basis for DFT models of <5-doped Si:P against which 
the localized basis techniques are assessed. Implications 
for the type of modeling being undertaken are discussed, 
and the models extended beyond those tractable with 
plane- wave techniques. Using these calculations, we ob- 
tain converged values for properties such as bandstruc- 
tures, energy levels, valley splitting, electronic densities 
of state and charge densities near the (5-doped layer. 

The paper is organized as follows: Sec. |TT] outlines 
the parameters used in our particular calculations; we 
present the results of our calculations in Sec. IIIII and 
conclusions are drawn in Sec. IIVI An elucidation of ef- 
fects modifying the bulk bandstructure follows in App. [A] 
&[B]to provide a clear contrast to the properties deriving 
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from the 5-doping of the silicon discussed in the paper. 



II. METHODOLOGY 

Density functional theory calculations have been car- 
ried out using both plane-wave and LCAO basis sets. 
For the plane- wave (PW) basis set, the Vienna ab ini- 
tio simulation package (vasp) 3 -!! software was used with 
projector augmented wave (PAW) 30 ' 31 pscudopotcntials 
for Si and P. Due to the nature of the plane-wave (PW) 
basis set, there exists a simple relationship between the 
cutoff energy and basis set completeness. For the struc- 
tures considered in this work, the calculations were found 
to be converged for PW cutoffs of 450 eV. 

Localized basis set calculations were performed us- 
ing the Spanish Initiative for Electronic Simulations 
with Thousands of Atoms (siesta)^ software. In this 
case, the P and Si ionic cores were represented by 
norm-conserving Troullier-Martins pseudopotentialsi^ 3 . 
The Kohn-Sham orbitals were expanded in the default 
single-C polarized (SZP) or double-^ polarized (DZP) ba- 
sis sets, which consist of 9 and 13 basis functions per atom 
respectively. Both the SZP and DZP sets contain s-, p-, 
and ci-type functions. These calculations were found to 
be converged for a mesh grid energy cutoff of 300 Ry. In 
all cases, the generalized gradient approximation (GGA) 
PBE 3 i exchange-correlation functional was used. 

The lattice parameter for bulk Si was calculated using 
an 8-atom cell, and found to be converged for all methods 
with a 12 x 12 x 12 Monkhorst-Pack (MP) fc-point mesh 3 ^. 
The resulting values are presented in Table U and were 
used in all subsequent calculations. 



Method ao 



PW (vasp) 5.469 
DZP (siesta) 5.495 
SZP (siesta) 5.580 



TABLE I. Equilibrium lattice parameters for an 8-atom cubic 
unit cell for the different methods used in this work. 

In modeling <5-doped Si:P, as used in Ref. [1(1 we 
adopted a tetragonal supercell description of the system, 
akin to that of Refs. & [lf| In accordance with ex- 
periment, we inserted the P layer in a monatomic (001) 
plane as one atom in four to achieve 25% doping. This 
will henceforth be referred to as 1/4 monolayer (ML) 
doping. In this case, the smallest repeating in-plane unit 
had four atoms per monolayer (to achieve 1 in 4 dop- 
ing), and was a square with sides parallel to the [110] 
and [110] directions. The square had side length a\[2 
(see Fig. [T]), where a is the simple cubic lattice constant 
of bulk silicon. The phosphorus layers had to be sepa- 
rated by a considerable amount of silicon due to the large 
Bohr radius of the hydrogen-like orbital introduced by P 
in Si (^2.5 nm). Ref. [l5| showed that this far exceeded 



the sub-nanometer cell side length. If desired, cells with 
a lower in-plane density of dopants may be constructed, 
by lengthening the cell in the x- and y-dircctions, such 
that more Si atoms occupy the doped monolayer in the 
cell. 




FIG. 1. (001) planar slice of the c (2 x 2) structure, at the 1/4 
ML doped monolayer. One of the Si sites has been replaced 
by a P atom (shown in dark gray). The periodic boundaries 
are shown in black. 

A collection of tetragonal cells comprised of 4, 8, 16, 32, 
40, 60, 80, 120, 160 and 200 monolayers were constructed, 
having four atomic sites per monolayer and oriented with 
faces in the [110], [110] and [001] directions (see Fig. 0. 
Cells used in PW calculations began at 4 layers and ran to 
80 layers; larger cells were not computationally tractable 
with this method. SZP and DZP models began at 40 
layers to overlap with PW for the converging region, and 
were then extended to their tractable limit (200 and 160 
layers, respectively) to study convergence past the capa- 
bility of PW. 




FIG. 2. Ball & stick model of a <5-doped Si:P layer, viewed 
along the [110] direction; 32 layers in the [001] direction are 
shown. Si atoms (small gray spheres), P atoms (large dark 
gray spheres), covalent bonds (gray sticks), repeating cell 
boundary (solid line). 

For the tetragonal cells the fc-point sampling was set 
as a 9 x 9 x JV T-centrcd MP mesh, as we have found that 
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failing to include T in the mesh can lead to anomalous 
placement of the Fermi level on bandstructure diagrams. 
N varied from 12 to 1 as the cells became more elongated 
(see Table HV] in App. [£}. 

Although it has been previously found that relaxing 
the positions of the nuclei gave negligible differences 
(< 0.005 A) to the geometry this was for a 12-layer 
cell and may not have included enough space between 
periodic repetitions of the doping plane for the full effect 
to be seen. We have performed a test relaxation on a 40- 
layer cell using the PW basis (vasp). The maximum sub- 
sequent ionic displacement was 0.05 A, with most being 
an order of magnitude smaller. The energy gained in re- 
laxing the cell was less than 37 meV (or 230 /xeV/atom). 
We therefore regarded any changes to the structure as 
negligibly small, and proceeded without ionic relaxation. 

Single-point energy calculations were carried out with 
both software programs; for vasp the electronic energy 
convergence criterion was set to 10 -6 eV, and the tetra- 
hedron method with Blochl correction^ was used. For 
SIESTA a two-stage process was carried out: Fermi-Dirac 
electronic smearing of 300 K was applied in order to con- 
verge the density matrix within a tolerance of 1 part in 
10 -4 ; the calculation was then restarted with smearing 
of K and a new electronic energy tolerance criterion of 
10" 6 eV applied (except for the 120- and 160-layer DZP 
models for which this was intractable; a tolerance of 10 -4 
eV was used in these cases). This two-stage process aided 
convergence as well as ensuring that the energy levels 
obtained were comparably accurate across methods. In 
addition, for each doped cell thus developed and studied, 
an undoped bulk Si cell of the same dimensions was con- 
structed to aid in isolating those features primarily due 
to the doping. 



III. RESULTS 
A. Analysis of bandstructure 

Once converged charge densities were obtained, band- 
structures were calculated along the M-T-X high- 
symmetry pathway (as shown in Fig. [TU]in Appendix [A"|) . 
using at least 20 Appoints between high-symmetry points. 
For comparative purposes, the bandstructures have all 
been aligned at the valence band maximum (VBM). 

Figure [3] contrasts bulk and doped bandstructures for 
the 40-layer PW calculation. DZP and SZP results are 
similar on this scale and are omitted in the interest of 
clarity in the diagram. As discussed in App. [B] it is ev- 
ident from the bulk values that the elongated cells have 
led to the folding of two CBM valleys towards the Y- 
point. Also visible is the difference the doping potential 
makes to the system; what was the lowest unoccupied 
orbital in the bulk is now dragged down in energy by the 
extra ionic potential. It is of note that the region near T, 
corresponding to the k z valleys which can be modeled as 
having different effective masses to the k x ^ y valleys j 1 ^ is 
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FIG. 3. Full band structure of the 40-layer tetragonal sys- 
tem, calculated using PW (vasp). Bulk and 1/4 ML doped 
structures are shown. 

brought lower than the region corresponding to the k XiV 
valleys and is non-degenerate. The second band behaves 
in a similar fashion. The third band appears to maintain 
a minimum away from the T-point in the SxET-direction 
(which is equivalent to the Apcc-direction; see App. |A"| . 
but in a less-parabolic fashion than the lower two; its 
minimum is similar to the value at T. This band is non- 
degenerate along this particular direction in £;-space, but 
due to the superccll symmetry it is actually 4-fold degen- 
erate, in contrast to the other bands. The Fermi level for 
the doped system is also shown, clearly being crossed by 
all three of these bands which are therefore able to act 
as open channels for conduction. 




0.25 M E T A 0.25 X 



FIG. 4. Bandstructure of the 40-layer tetragonal system, 
zoomed in on the A band, calculated with (a) PW (vasp), 
(b) DZP (siesta), and (c) SZP basis sets. 

As mentioned above, the bandstructures are similar 
across all methods, but upon detailed inspection impor- 
tant differences come to light. A closer look at the A 
band shows a qualitative difference between the prcdic- 
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tions using SZP (Fig. gfc) and the PW and DZP re- 
sults (Figs. &Hb) : the models with a more complete 
basis predict the band minimum to occur in the Etet 
(Apcc) direction, below the value at T, while the SZP 
bandstructure shows the reverse - the minimum at T, a 
similar amount below a secondary minimum in the £tet 
direction. 
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The difference between the energies of the first two 
band minima (Fi-T2, illustrated in Fig. [5]), or the valley 
splitting, from the PW and DZP calculations agree with 
each other to within ~6 meV. Significantly, the value ob- 
tained using the SZP basis set differs by 52 meV, some 
55% larger than the value obtained using the PW ba- 
sis set. The importance of this discrepancy cannot be 
overstated; this valley splitting is directly relatable to 
experimentally observable resonances in transport spec- 
troscopy of devices made with this (5-doping technology 
(see Ref. 

In the smallest cells (< 16 layers), less than three bands 
are observed. This is likely due to the lack of cladding 
in the z-direction, leading to significant interaction be- 
tween the dopant layers, raising the energy of each band. 
Whilst the absolute energy of each level still varies some- 
what, even with over 100 layers incorporated, we find 
that the ri-T^ values are well- converged with 80 layers 
of cladding for all methods (see Fig. [5|). Indeed, they may 
be considered reasonably converged even at the 40-layer 
level (0.5 meV or less difference to the largest models 
considered). The differences between the energies of the 
second and third band minima (r 2 -A splittings) are also 
shown in Fig. [5J and show good convergence (within 1 
meV) for cells of 80 layers or larger. 

The Fermi level follows a similar pattern to the T- and 
A-levels. In particular, the gap between the Fermi level 
and Ti level does not change by more than 1 meV from 
60 to 160 layers. 

Given that the properties of interest are the differences 
between the energy levels, rather than their absolute val- 
ues (or position relative to the valence band), in the in- 
terest of computational efficiency we observe that using 
the DZP basis with 80 layers of cladding is sufficient to 
achieve consistent, converged results. 



B. Valley splitting 
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FIG. 5. Minimum band energies for tetragonal systems with 
1/4 ML doping, calculated using : (a) PW (vasp), (b) DZP 
(siesta), and (c) SZP (siesta) basis sets. Fermi level also 
shown where appropriate. Bold numbers indicate energy dif- 
ferences between band minima. 



Tabic pi] summarizes the valley splitting values of 1/4 
ML P-dopcd silicon obtained using different techniques, 
showing a large variation in the actual values. In order 
to make sense of these results, it is important to note 
two major factors that affect valley splitting: the doping 
method and the arrangement of phosphorus atoms in the 
(5-layer. As the results from Ref. [l6| show, the use of im- 
plicit doping causes the valley splitting value to be much 
smaller than in an explicit case (~7 meV vs. 120 meV). It 
is also shown that the use of random P coverage on the 5- 
layer reduces the valley splitting value by only 40-50 meV 
compared to the fully ordered placement, leaving a large 
discrepancy between the valley splitting results from im- 
plicit and explicit doping. This massive decrease in valley 
splitting due to implicit doping can be explained by the 
smearing of the doping layer in the direction normal to 
the (5-layer, thereby decreasing the quantum-confinement 
effect responsible for breaking the degeneracy in the sys- 
tem. Ref. [l6| also shows that the arrangement of the 



phosphorus atoms in the 5-layer strongly influences the 
valley splitting value. In particular, they showed that 
there is a difference of up to 220 meV between P doping 
along the [110] direction and along the [100] direction. It 
must be noted, however, that most of their patterns are 
not yet physically realizable due to the P incorporation 
mechanism currently employed. 
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a implicit doping 
b explicit doping 
c M X M X 1 fc-points used 

d M X M X N fc-points used; TV as per Table llVl in App. |X] 

TABLE II. Valley splitting values of 1/4 ML P-doped silicon 
obtained using different techniques. Techniques are grouped 
by similarity. 

Our results show that valley splitting is highly sensi- 
tive to the choice of basis set. Due to the nature of PW 
basis set, it is straightforward to improve its complete- 
ness by increasing the plane wave cutoff energy. In this 
way, we establish the most accurate valley splitting value 
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within the context of density functional theory. Using 
this benchmark value, we can then establish the validity 
and accuracy of other basis sets, which can be used to 
extend the system sizes to that beyond what is practical 
using PW basis set. As seen in Table [Til the valley split- 
ting value converges to 93 meV using 80-layer cladding. 
The DZP localized basis set gives an excellent agreement 
at 99.5 meV using 80-laycr cladding (representing a 7% 
difference). On the other hand, the SZP localized basis 
set (similar to what was used in Refs. [H and HH) gave 
a value of 145 meV using the same amount of cladding. 
This represents a significant difference of 55% over the 
value obtained using PW basis set, and demonstrates 
that the SZP basis set is unsuitable for accurate deter- 
mination of valley splitting in these systems. 



C. Density of states 

The electronic density of states (cDOS) was calculated 
for each cell. Figure |5] compares the unsealed cDOS for 
bulk 80-layer cells to that of doped cells varying from 
40 to 80 layers. The bulk bandgap is visible, with the 
conduction band rising sharply to the right of the figure. 
The doped eDOS exhibits density in the bulk bandgap, 
although the features of the spectra differ slightly accord- 
ing to the basis set used. 




FIG. 6. Electronic densities of states for tetragonal systems 
with and 1 /4 ML doping, calculated using the DZP (siesta) 
basis set. The Fermi level is indicated by a solid vertical line 
with label. 50 meV smearing was applied for visualization 
purposes. 

The Fermi energy exhibits convergence with respect 
to the amount of cladding, as reported above. It is also 
notable that the eDOS within the bandgap are nearly 
identical regardless of the cell length (in z). This indi- 
cates that layer-layer interactions are negligibly affect- 
ing the occupied states, and therefore that the applied 
"cladding" is sufficient to insulate against these effects. 
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D. Electronic width of the plane 

In order to quantify the extent of the donor-electron 
distribution, we have integrated the local density of states 
(LDOS) between the VBM and Fermi level and taken the 
planar average with respect to the z-position. Figure [7] 
shows the planar average of the donor electrons (a sum 
of both spin-up and spin-down channels) for the 80-layer 
cell calculated using the DZP basis set. After removing 
the small oscillations related to the crystal lattice to focus 
on the physics of the 5-layer, by Fourier transforming, a 
Lorentzian function was fitted to the distribution profile. 
(Initially, a three-parameter Gaussian fit similar to that 
used in Ref. [24| was tested, but the Lorentzian gave a 
better fit to the curve.) 
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FIG. 7. Planar average of the donor-electron density as a 
function of z-position for the 1/4 ML doped, 80-layer cell 
calculated using the DZP basis set. The fitted Lorentzian 
function is also shown. 

Table IIIII summarizes the maximum donor-electron 
density and the full width at half maximum (FWHM) 
for the 1/4 ML doped cells, each calculated from the 
Lorentzian fit. Both of these properties are remarkably 
consistent with respect to the number of layers, indicat- 
ing that they have converged sufficiently even at 40 lay- 
ers. 
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layers 


(xlO" 3 e/A) 
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40 


3.8 


6.2 


60 


3.9 


6.2 


80 


3.9 


6.5 



TABLE III. Calculated maximum donor-electron density, 
Pmax, and full width at half maximum, FWHM, as a func- 
tion of the number of layers in the 1/4 ML doped cells. The 
DZP basis was used. 

Our results differ from a previous DFT calculation-^ 
which cited a FWHM of 5.6 A for a 1/4 ML doped, 
80-layer cell calculated using the SZP basis set (and 
10x10x1 /c-points). We note that those values were 
taken from the unfitted, untransformed donor-electron 



distribution, and represent a ^15% underestimation of 
the DZP result. 



IV. CONCLUSION 

In this article, we have studied the valley splitting of 
monolayer <5-doped Si:P, using a density functional theory 
model with a plane-wave basis to establish firm grounds 
for comparison with less computationally-intensive local- 
ized basis ab initio methods. We found that the best cur- 
rent descriptions of these systems (by density functional 
theory, using SZP basis functions) overestimate the val- 
ley splitting by over 50%, due to an assumption made 
early in their methodology. We show that DZP basis 
sets are complete enough to deliver values within 10% of 
the plane-wave values, and due to their localized nature, 
are capable of calculating the properties of models twice 
as large as is tractable with plane-wave methods. These 
DZP models are converged with respect to size well be- 
fore their tractable limit, which approaches that of SZP 
models. 

Valley splittings are important in interpreting trans- 
port spectroscopy experiment data, where they relate to 
families of resonances, and in benchmarking other theo- 
retical techniques more capable of actual device model- 
ing. It is therefore pleasing to have an ab initio descrip- 
tion of this effect which is fully-converged with respect to 
basis completeness, as well as the usual size effects and 
fc-point mesh density. 

We have also studied the bandstructures with all three 
methods, finding that the DZP correctly determines the 
A-band minima away from the T point, where the SZP 
method does not. We show that these minima occur in 
the £ direction for the type of cell considered, not the 
A direction as has been previously reported. Having es- 
tablished the DZP methodology as sufficient to describe 
the physics of these systems, we then calculated the elec- 
tronic density of states, and the electronic width of the 5- 
layer. We found that previous SZP descriptions of these 
layers underestimate the width of the layers by almost 
15%. 

We have shown that the properties of interest of i5- 
doped Si:P are well-converged for 40-layer supercclls us- 
ing a DZP description of the electronic density. We rec- 
ommend the use of this amount of surrounding silicon, 
and technique, in any future DFT studies of these and 
similar systems - especially if inter-layer interactions are 
to be minimized. 
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Appendix A: Subtleties of bandstructure 

Regardless of the type of calculation being undertaken, 
a bandstructure diagram is inherently linked to the type 
(shape and size) of cell being used to represent the sys- 
tem under consideration. For each of the 14 Bravais lat- 
tices available for three-dimensional supercclls, a particu- 
lar Brillouin zone (BZ) with its own set of high-symmetry 
points exists in reciprocal spaced. Similarly, each BZ has 
its own set of high-symmetry directions. Some of these 
BZs share a few high-symmetry point labels (or direc- 
tions), such as X or L (A or E), and they all contain T, 
but these points are not always located in the same place 
in reciprocal space. 

A simple effect of this can be seen by increasing the 
size of a supercell. This has the result of shrinking the 
BZ, and the coordinates of high-symmetry points on its 
boundary, by a corresponding factor. Consider the con- 
duction band minimum (CBM) found at the A valley 
in the Si conduction band. This is commonly located 
at kn ~ 0.85— in the A direction towards X. Should 

u a 

we increase the cell by a factor of 2, the BZ will shrink 
(BZ— >BZ ), placing the valley outside the new BZ bound- 
ary (past X ); but a valid solution in any BZ must be 
a solution in all BZs. This results in the phenomenon 
of band folding, whereby a band continuing past a BZ 
boundary reenters the BZ on the opposite side. Since the 
X direction in a face-centred cubic (FCC) BZ is 6-fold 
symmetric, a solution near the opposite BZ boundary is 
also a solution near the one we are focussing on. This 
results in the appearance that the band continuing past 
the BZ boundary is "reflected", or folded, back on it- 
self into the first BZ. Since the new BZ boundary in this 
direction is now at fc BZ = X =0.5—, the location of 

the valley will be at k' = X - (k - X^j ~ 0.15^, as 

mentioned in Ref. [IH Each further increase in the size 
of the supercell will result in more folding (and a denser 
bandstructure). Care is therefore required to distinguish 
between a new band and one which has been folded due 
to this effect when interpreting bandstructure. 

Continuing with our example of silicon, whilst the clas- 
sic bandstructure^ is derived from the bulk Si primitive 
FCC cell (containing two atoms), it is often more con- 
venient to use a simple cubic (SC) supercell (8 atoms) 
aligned with the (100) crystallographic directions. In this 
case, we experience some of the common labelling; the A 
direction is defined in the same manner for both BZs, al- 
though we see band folding (in a similar manner to that 
discussed above) due to the size difference of the recipro- 
cal cells (see Fig. [5]). We also see a difference in that al- 
though the E direction is consistent, the points at the BZ 
boundaries have different symmetries and therefore, la- 
bels (if fcc , Mac)- (The £FCC-P om t and Apcc-direction 
have no equivalent for tetragonal cells, and hence we do 
not consider bandstructure in that direction here) 

Consider now the (5-doping case discussed above (see 
Sec. HI]) , where we wish to align our cell with the [110] 



and [110] directions (by rotating the cell 45° anticlock- 
wise about z; this will also require a resizing of the cell 
in the plane to maintain periodicity - see Fig. [5]), to al- 
low us to include precisely four atoms per monolayer (as 
required for the minimal representation of 1/4 ML dop- 
ing). We now have a situation where the Xtet point in 
the new tetragonal BZ (see Fig. fTUj) is no longer in the 
direction of the X^c point in the simple cubic BZ, de- 
spite both X points being in the centre of a face of their 
BZ. Due to the rotation, what was the Age direction in 
the simple cubic BZ is now the Etet direction (pointing 
towards M, at the corner of the BZ in the k z = plane) 
in the tetragonal BZ. The tetragonal CBM, while phys- 
ically still the same as the CBM in the FCC or simple 
cubic BZ, is not represented in the same fashion (see Fig. 
EU). 



Appendix B: Band folding in z-direction 

Increasing the z-dimension of the cell leads to succes- 
sive folding points being introduced as the Brillouin zone 



=0 

a 





(b) 



(c) 



FIG. 8. (a) Typical band structure of bulk Si for 2-atom FCC 
(solid lines) and 8-atom SC cells (dotted lines with squares), 
calculated using the VASP plane- wave method (see Sec. HTjl . 
(b) 2-atom FCC cell, (c) 8-atom SC cell. 
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c(2x2) 



p(lxl) 



p(2x2) 



FIG. 9. Geometrical difference between the simple cubic and 
tetragonal cells; (001) planar cut through an atomic mono- 
layer. 




FIG. 10. The Brillouin zone for a tetragonal cell. The M-T- 
X path used in this work is shown. 




FIG. 11. Band structure diagram for the tetragonal bulk Si 
structures with increasing number of layers, calculated using 
the vasp plane wave method (see Sec. [TTJ . 



(BZ) shrinks along k z (see App. IA|) . This has the effect 
of shifting the conduction band minima in the ±fc z direc- 
tions closer and closer to the T point (see Fig. 8(a) ) and 



making the bandstructurc extremely dense when plotting 
along k z . This results in the value of the lowest unoccu- 
pied eigenstate at T being lowered as what were originally 
other sections of the band are successively mapped onto 
r, and after a sufficient number of folds the value at T 
is indistinct from the original conduction band minimum 
(CBM) value. The effects of this can be seen in Table HVl 



which describes increasingly elongated tetragonal cells of 
bulk Si. When we then plot the bandstructurc in a differ- 
ent direction, e.g. along k x , the translation of the minima 
from ±fc 2 onto the F-point appear as a new band with 
two- fold degeneracy. The degeneracy of the original band 
drops from 6- to 4-fold, in line with the reduced symme- 
try (we only explicitly calculate one, and the other three 
occur due to symmetry considerations) . This is the origin 
of the 'T-bands" discussed in Refs. [lj &[l5|. Once the k z 
valleys are sited at T, parabolic dispersion corresponding 
to the transverse kinetic energy terms is observed along 
k x and k y , at least close to the band minimum (see Fig. 
□I]). 

All methods considered in Table |IV] show the LUMO 
at r (folded in along ±k z ) approaching the CBM value 
as the amount of cladding increases; at 80 layers, the 
LUMO at T is within 1 mcV of the CBM value. It is also 
of note that the PW indirect bandgap agrees well with 
the DZP value, and less so with the SZP model. This is 
an indication that, although the behaviour of the LUMO 
with respect to the cell shape is well-replicated, the SZP 
basis set is demonstrably incomplete. Conversely, pair- 
wise comparisons between the PW and DZP results show 
agreement to within 5 meV. 

It is important to distinguish effects indicating con- 
vergence with respect to cladding for doped cells (i.e. 
elimination of layer-layer interactions) from those men- 
tioned above which derive from the shape and size of 
the supercell. Strictly, the convergence (with respect to 
the amount of encapsulating Si) of those results we wish 
to study in detail, such as the differences in energy be- 
tween occupied levels in what was the bulk bandgap, pro- 
vide the most appropriate measure of whether sufficient 
cladding has been applied. 



Basis 


No. of No. of LUMO 


CBM 


type 


layers 


fc-pts 


at r 


(at Afcc) 






in k z 


(eV) 


(eV) 


PW 


4 


12 


0.7517 




(VASP) 


8 


6 


0.7517 






16 


3 


0.6506 






32 


2 


0.6170 






40 


1 


0.6179 






64 


1 


0.6137 






80 


1 


0.6107 


0.6102 


DZP 


40 


1 


0.6218 




(siesta) 


60 


1 


0.6194 






80 


1 


0.6154 






120 


1 


0.6145 






160 


1 


0.6151 


0.6145 


SZP 


40 


1 


0.8392 




(siesta) 


60 


1 


0.8349 






80 


1 


0.8315 






120 


1 


0.8311 






160 


1 


0.8315 






200 


1 


0.8310 


0.8309 



TABLE IV. Energy levels of tetragonal bulk Si structures. 
(For details of calculation parameters, see Sec. [TTJ) 
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